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Marine  biological  models  are  usually  complex  with  many  free  parameters.  Parameter  pri¬ 
oritization  (based  on  contribution  to  model  output)  is  important  for  system  management 
but  difficult.  A  variance-based  sensitivity  analysis  is  developed  in  this  paper  using  the 
Sobol’-Saltelli  sensitivity  indices,  which  measure  the  relative  importance  of  each  parameter 
(or  group  of  parameters)  and  range  these  parameters  along  their  contribution  to  output  vari¬ 
ability.  To  reduce  the  number  of  degrees  of  freedom,  the  model  output  is  decomposed  using 
the  warping  functions  or  irreversible  predictability  time.  A  simple  three-component  [nutri¬ 
ents,  phytoplankton  and  zooplankton  (NPZ)]  model  with  23  parameters  for  reproducing 
annual  phytoplankton  cycle  of  the  Black  Sea  is  taken  as  the  example  to  show  the  usefulness 
and  procedure  of  the  sensitivity  analysis.  Single  and  total  sensitivity  indices  showed  strong 
sensitivity  of  the  biological  model  to  the  light  limitation  of  the  phytoplankton  growth.  This 
agrees  well  with  physical  intuition.  However,  ranging  model  parameters  along  their  con¬ 
tributions  to  model  output  variability  does  not  follow  exactly  the  physical  intuition  when 
model-related  errors  from  large  perturbations  of  the  parameters  are  not  small.  For  example, 
the  model  output  becomes  very  sensitive  to  the  nutrient  stock  parameterization  for  certain 
combinations  of  the  light-related  factors. 

©  2007  Elsevier  B.V.  All  rights  reserved. 


1.  Introduction 

Ocean  models,  especially  ocean  biological  models,  in  general, 
have  many  uncertain  parameters,  which  should  be  identi¬ 
fied  from  data  or  the  physics  (Lozano  et  al.,  1996;  Omlin 
et  al.,  2001;  Fulton  et  al.,  2004;  Lermusiaux  et  al.,  2006 
among  others).  Various  data  assimilation  methods  may  be 
used  for  model  parameter  identification:  the  adjoint  method 
(Evensen  et  al.,  1998),  the  non-linear  optimization  technique 
(Fasham  and  Evans,  1995),  the  weak-constraint  parameter 
estimation  (Loza  et  al.,  2004)  and  others.  The  basic  con¬ 
cept  of  these  methods  is  to  vary  model  parameters  until 
the  misfit  between  temporally  varying  modeled  and  observed 


data  is  minimal,  while  the  model  equations  are  satisfied 
exactly. 

Although  robust  dynamical  regimes  (attractors)  repro¬ 
duced  by  biological  models  are  not  very  complex  (most  of 
such  models  demonstrate  only  simple  periodical  or  quasi¬ 
periodical  behavior),  the  parameter  identification  is  quite  a 
difficult  problem  by  a  number  of  reasons.  First,  data  and  model 
may  be  incompatible  because  the  data  contain  contributions 
from  hydrodynamic  and  biological  processes,  which  may  not 
be  resolved  by  the  model.  Model  error  (no  matter  how  small  it 
is)  can  cause  the  solutions  deviating  far  from  the  data.  For 
example,  Fasham  and  Evans  (1995)  could  not  find  a  single 
parameter  set  that  fits  the  observational  data  well.  Spitz  et 
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al.  (1998)  could  not  estimate  the  optimal  model  parameters 
using  observational  data. 

Second,  modern  data  assimilation  methods  are  on  the  base 
of  statistical  estimation  theory,  which  uses  the  same  founda¬ 
tion  as  the  Kalman  filtering:  (1)  data  and  model  are  assumed 
to  be  unbiased;  (2)  error  variances  and  co-variances  (dictating 
model-data  difference)  are  used  to  correct  the  model  state; 
(3)  Gaussian  statistics  are  assumed  for  the  data  errors.  These 
conditions  may  not  be  true  in  biological  modeling. 

Third,  even  simple  biological  models  with  3-10  model  vari¬ 
ables  often  contain  20-30  or  more  model  parameters  (e.g., 
Fasham  et  al.,  1990;  Oguz  et  al.,  1996;  Fennel  et  al.,  2001; 
Kantha,  2004).  When  the  number  of  model  variables  is  con¬ 
siderably  less  than  the  number  of  model  parameters,  the 
parameter  identification  does  not  have  a  unique  solution, 
and  in  general,  selection  of  the  “optimal”  solution  is  difficult. 
This  leads  to  the  identifiability  problem  in  ecological  model¬ 
ing,  which  concerns  the  uniqueness  of  the  model  parameters 
determined  from  the  input-output  data,  under  ideal  condi¬ 
tions  of  noise-free  observations  and  error- free  model  structure 
(Beck,  2002).  The  statistical  method  identifies  the  model 
parameters  with  only  minimum  difference  between  model 
and  data  but  does  not  guarantee  the  absolute  minimum  error 
(Schittkowski,  2002):  hyper-parameterized  models  may  have 
many  possible  solutions. 

In  many  cases,  change  in  certain  parameters  (non-control 
parameters)  causes  only  little  change  in  model  output.  There¬ 
fore,  these  parameters  can  be  approximately  determined  and 
then  fixed.  Change  in  other  parameters  (control  parameters) 
causes  large  change  in  model  output.  Thus,  the  control  param¬ 
eters  have  to  be  determined  in  a  very  accurate  manner  because 
they  affect  the  model  predictability.  The  question  arises,  how 
can  we  range  model  parameters  according  to  their  contribu¬ 
tions  to  model  output  variability? 

Here,  a  phenomenological  approach  may  be  used  to 
detect  model  sensitivity  to  the  control  parameter[s].  Such  an 
approach  requires  rich  practical  experience  and,  in  general, 
often  gives  reasonable  results  if  the  number  of  control  param¬ 
eters  is  not  large.  If  the  number  of  model  parameters  is  large 
(marine  biological  models  are  the  case),  such  an  approach  may 
fail  no  matter  how  rich  a  researcher’s  experience  is,  because 
model  sensitivity  relative  to  one  parameter  often  differs  from 
model  sensitivity  to  a  group  of  parameters. 

Alternative  methods  to  determine  the  control  parameters 
are  the  first-order  sensitivity  function  (Chu  et  al.,  2004)  and 
the  adjoint  method  (Evensen  et  al.,  1998).  The  traditional  sen¬ 
sitive  analysis  based  on  the  direct-perturbation  method  (e.g., 
Dickinson  and  Gelinas,  1976)  is  popular  in  biological  oceanog¬ 
raphy.  For  example,  Oguz  et  al.  (1996)  used  this  approach  to 
verify  a  low-component  model  of  annual  phytoplankton  cycle 
in  the  Black  Sea.  The  direct-perturbation  method  pursuant 
to  which  the  sensitivity  of  model  output  to  change  in  model 
parameters]  is  found  by  comparing  model  integrations  with 
the  only  (finite)  difference  in  the  parameter  of  interest.  The 
disadvantage  of  the  direct  method  is  that  separate  model  inte¬ 
gration  must  be  performed  for  each  parameter  of  interest. 
That  a  priori  assumes  additive  contribution  from  each  param¬ 
eter  to  model  output. 

The  adjoint  method  (Lawson  et  al.,  1995;  Evensen  et  al., 
1998  and  others)  estimates  model  parameters  and  variables 


through  fitting  the  model  to  data,  using  model  equation  as  a 
constraint.  However,  the  method  requires  an  initial  guess  for 
unknown  initial  conditions  and  parameters.  Second,  a  biolog¬ 
ical  model  cannot  be  taken  as  a  ‘true  model’  because  of  many 
parameterization  schemes  involved.  Third,  although  Lawson 
et  al.  (1995)  reported  that  the  adjoint  method  worked  reason¬ 
ably  well  even  for  “data”  with  20%  noise-to-signal  level,  it  is  not 
clear  how  the  optimal  model  parameters  are  determined.  Pires 
et  al.  (1996)  pointed  out  that  for  non-linear  dynamical  mod¬ 
els  and  noisy  data  there  are  limitations  in  application  of  the 
adjoint  technique,  and  its  convergence  to  the  optimal  solution 
is  not  obvious  with  the  presence  of  noise  in  the  data. 

As  the  exact  values  of  control  parameters  of  a  biolog¬ 
ical  model  are  unknown,  the  linear  sensitivity  approach 
assumes  explicitly  no  interactions  among  forecast  model- 
related  errors  caused  by  parameter  perturbations.  In  many 
practical  cases,  this  assumption  is  unrealistic,  and  the  model 
regimes  and  transitions  among  regimes  are  controlled  by 
parameters  determined  from  the  sensitivity  analysis  on  finite- 
amplitude  parameter  perturbations  (Nicolis,  2003). 

The  primary  goals  of  the  proposed  study  are  outlined  as  fol¬ 
lows:  (a)  develop  a  model-independent  non-linear  sensitivity 
analysis  for  marine  biological  models  using  the  Sobol’-Saltelli 
sensitivity  indices  (Saltelli  et  al.,  1993,  2000,  2005).  (b)  Use  spe¬ 
cial  metrics,  such  as  warping  functions  and  the  irreversible 
predictability  time  (IPT)  (Chu  et  al.,  2002a, b,c)  as  model  out¬ 
put.  IPT  is  developed  on  the  base  of  first  passage  time,  (c) 
Demonstrate  capability  of  this  approach  through  the  analysis 
of  a  three-component  (nutrients,  phytoplankton,  zooplank¬ 
ton)  model  for  the  annual  phytoplankton  cycle  in  the  Black 
Sea.  The  choice  of  the  model  is  from  research  interests  of  the 
authors,  and  is  not  principle. 

The  non-linear  sensitivity  analysis  does  not  find  the  opti¬ 
mal  model  parameters  directly.  It  assesses  the  influences  or 
relative  importance  of  each  model  parameter  to  the  model 
output  and  determines  which  parameters  are  control  param¬ 
eters  contributing  most  to  the  output  variability  and,  possibly, 
requiring  additional  research  to  reduce  output  uncertainty, 
and  which  parameters  are  non-control  parameters  and  can  be 
estimated  approximately.  Excluding  the  non-control  parame¬ 
ters,  we  may  reduce  the  number  of  model  parameters  that  are 
identified  from  data  or  physics. 

The  rest  of  the  paper  is  organized  as  follows.  Section 
2  describes  the  sensitivity  analysis  using  the  Sobol’-Saltelli 
indices.  Section  3  presents  model  output  representations  for 
estimating  the  Sobol’-Saltelli  indices.  Section  4  shows  the 
model  output  representation  using  the  warping  functions  and 
IPT  for  the  non-linear  sensitivity  analysis.  Section  5  describes 
the  simplified  three-component  biological  model  for  the  Black 
Sea  phytoplankton  annual  cycle  (hereafter,  the  NPZ  model). 
Section  6  depicts  the  experiment  design.  Sections  7-10  present 
the  results  and  their  oceanographic  interpretations.  Section  11 
presents  the  conclusions. 


2.  Non-linear  sensitivity 

A  variance-based  method  (Saltelli  et  al.,  2000,  2005)  is  devel¬ 
oped  to  estimate  the  non-linear  sensitivity  of  a  biological 
model  to  large  variations  of  model  parameters.  Following 
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Sobol’  (1993),  Saltelli  et  al.  (1993)  and  Saltelli  (2002),  we  intro¬ 
duce  the  sensitivity  indices,  which  are  quantitative  measures 
of  model  output  sensitivity  to  large  variations  of  model  param¬ 
eters. 

For  simplicity,  consider  a  scalar  model  output  y=/(ai,  a2, 
. . .,  aje)  corresponding  to  a  number  of  non-correlated  model 
parameters  a  =  (ai ,  a2 , . . . ,  a^)  with  the  j  oint  probability  density 
function: 

k 

P(ai,  a2 . ak)  =  J]pi(a.).  (2.1) 

i=l 

Here,  p; (a;)  is  the  PDF  for  ith  parameter. 

The  total  sensitivity  of  the  model  output  relative  to  vari¬ 
ations  of  all  model  parameters  is  estimated  by  the  variance 
V 

k 

[f(«i.  “2,  ■  ■  ■ ,  ak)  -  E(y)]2  J]pi(aj)  da; 

i=l 


where  a _j  =  (&i,  a2 _ _  aj_1:  aJ+1, _ afe)  is  the  sub-vector  of  a 

containing  all  the  varying  parameters  other  than  ap 

E(y|a_j)  = 

1  f  &  1,  a2, ....  ap  ... ,  ak)pj(Uj )  daj. 

(2.9) 

u-;(y)  =  j 

n  fe 

E2(yla-j)  Pifii)  dfij. 

i  =  1 

¥i 

(2.10) 

The  total  sensitivity  of  the  model  output  to  variations  of 
parameter  Oj  is  written  by 

V-VT 

S j=^-L,  (2.11) 

which  shows  the  non-additive  part  of  model  output  sensitiv¬ 
ity  caused  by  interactions  among  model-related  errors  due  to 
perturbations  of  model  parameters.  Following  Sobol’  (1993), 
the  total  output  variance  V  is  decomposed  by: 


=  E(y2)  —  E2(y),  (2.2) 

where  E(y)  =  //.../  f(a1:  a2, . . . ,  afe)nLiPi(ai)  dai- 

The  output  variance  for  the  fixed  model  parameter  (qj), 
should  be  determined  by  the  conditional  variance  as 

Vj  =  V(y)  -  E{V(y|Oj)}  =  V{E(y|flj)}  =  Uj  -  E2(y),  (2.3) 

where  E(y|aj)  is  the  conditional  mean: 


E(y|flj) 


f(ai,a2,. 


■■ak)  ^Qpi(ai)  daj. 


i  =  1 

¥j 


(2.4) 


Uj  = 


E2(yl aj)Pj(aj)  da,-. 


(2.5) 


The  conditional  variance  Vj  is  a  good  measure  of  the  sen¬ 
sitivity  of  y  with  respect  to  the  parameter  (aj).  Once  divided  by 
the  unconditional  variance  V,  it  is  called  first-order  (or  single) 
sensitivity  indices 


Sj  = 


v(y)' 


(2.6) 


fe 

V  =  XVJ  +  XVjJ  +  XVi>'  +  ■  ■  ■  +  vi2...fe-  (2-12) 

j=l  j<i  j<i<l 

which  shows  partition  of  the  variance  between  the  main 
effects  (defined  by  Vj)  and  the  interaction  terms  (defined  by 
Vy>  Vs,,  . . .).  It  is  clear  that  when  Eq.  (2.12)  holds,  we  can  iden¬ 
tify  V_j  =  V{(y|a  j)}  as  the  sum  of  all  terms  in  the  right-hand 
side  of  (2.12)  without  terms  including  the  subscript  j. 

Therefore,  for  the  total  sensitivity  case,  the  index  for  the 
first  parameter  of  a  three-parameter  model  is 

Sj  =  Si  +  Su  +  Si3  +  S123.  (2.13) 

Dividing  both  parts  of  (2.12)  by  V,  we  get 
fe 

5>  +  _|_  £5*  +  ■  ■  ■  +  Si2...k  =  1,  (2.14) 

j= 1  j<i  j<  id 

which  shows  that  any  sensitivity  index  Sj  varies  between  0 
and  1. 

Sj  =  1 

for  the  model  output  depending  only  on  the  parameter  cij,  and 


Similarly  to  (2.6)  second-order  sensitivity  indices  deter¬ 
mining  model  output  sensitivity  relative  to  variations  of  two 
parameters  as  and  oij  (j  >  i )  may  be  introduced  by 

Sij  =  t7^)  •  V'J  =  V{E(ylQD  Q;))  -  v>  -  vi .  (2-7) 

which  show  the  joint  effects  of  both  parameters  af  and  a.j  on 
the  model  output. 

The  model  sensitivity  relative  to  variations  of  all  the 
parameters  excluding  a.j,  is  represented  by  the  conditional 
variance, 

(2.8) 


Sj  =  0,  sj  =  0 

when  the  model  output  does  not  depend  on  the  parameter  ap 
It  is  also  clear  from  Eqs.  (2.11)  and  (2.14)  that 

XSI  5  1  5  XSJ’  <2-15> 

1=1  J 

with  equalities  only  when  all  interaction  terms  in  (2.14)  are 
zero. 

The  major  advantage  of  variance-based  non-linear  analy¬ 
sis  is  that  we  can  account  for  and  specify  the  contributions 
of  model-related  errors  caused  by  perturbations  of  differ¬ 
ent  parameters  and  interactions  among  these  errors  using 


Vj  =  V{E(y | a  p)}  =  U_j(y)  -  E2(y), 
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second-order  indices  S;;,  Sj,i,  ....  In  this  case,  the  total  model 
output  is  a  non-additive  function  of  model  parameters. 

Use  of  the  sensitivity  indices  requires:  (a)  specifying  the 
model  output  y  in  the  simplest  form  and  (b)  developing  a 
numerical  technique  to  estimate  appropriate  variances  of  the 
model  output.  If  the  Monte-Carlo  method  (the  approach  used 
in  the  present  study)  is  used  to  compute  the  conditional  means 
and  variances,  the  model  output  should  be  represented  in  the 
simplest  form  so  that  the  sensitivity  analysis  would  be  com¬ 
putationally  feasible. 


3.  Model  output  representation 

We  suggest  to  use  the  following  representations  for  model 
output.  First,  considering  the  non-negative  feature  of  output 
for  biological  ocean  models,  the  warping  functions  (Gervini 
and  Gasser,  2004)  can  be  used  as  for  identification  of  func¬ 
tion  landmarks,  such  as  the  phytoplankton  blooms.  Second, 
IPT  introduced  by  Ivanov  et  al.  (1994),  Chu  et  al.  (2002a, b,c) 
and  Chu  and  Ivanov  (2005)  to  quantify  the  model  predictabil¬ 
ity  for  finite-amplitude  errors,  is  a  natural  measure  of  model 
sensitivity  when  perturbations  of  model  parameters  are  not 
small. 

3.1.  Warping  functions 

Following  Gervini  and  Gasser  (2004),  a  non-negative  one- 
variable  model  output  (y(t)>0,  such  as  phytoplankton 
concentration)  is  represented  by 

y(t)  =  sF{u(t)}  +  5(t),  (3.1) 


The  parameters  of  (3.1)— (3.3)  are  estimated  by  minimizing 
the  mean  integrated  squared  error, 

rt o+T  \ 

Jt  Iy(t)  -  sFju(t)}]2  dty  -»  min,  (3.4) 

where  ()  is  averaging  over  a  statistical  ensemble  of  the 
model  output.  Gervini  and  Gasser  (2004)  suggest  a  two- 
stage  algorithm  (a  MATLAB  routine  is  available  on  web 
page  http://www.unizh.ch/biostat/People/gervini)  to  mini¬ 
mize  (3.4),  and  point  out  that  only  a  few  functions  {Vq(t)}  are 
often  needed  to  reconstruct  geometrical  specificities  of  the 
model  output,  y(t)  >  0.  Therefore,  only  few  parameters  s  and  dj 
can  be  used  to  estimate  non-linear  sensitivity  model  output 
for  a  biological  model. 

3.2.  IPT 

The  IPT  is  defined  as  a  time-period  when  the  prediction  error 
first  exceeds  a  pre-determined  criterion  (i.e.,  the  tolerance 
level).  The  probability  density  function  of  IPT  with  a  given 
initial  error  satisfies  the  backward  Pontryagin-Kolmogorov- 
Stratonovich  equation.  Using  IPT  as  a  quantitative  measure  for 
prediction  skill,  both  linear  and  non-linear  regimes  of  forecast 
errors  were  found  in  the  low-order  atmospheric  model  (Chu 
et  al.,  2002a, b)  and  regional  ocean  circulation  model  (Chu  et 
al.,  2002c;  Chu  and  Ivanov,  2005). 

Following  Ivanov  et  al.  (1994)  IPT  (r)  is  determined  as  a  time 
for  which  difference  between  individual  model  output  z  and 
some  reference  solution  z  (defined  as  a  model  solution  with 
specified  model  parameters  and  initial  conditions)  will  exceed 
the  tolerance  level  (or  the  accepted  prediction  error)  s  for  the 
first  time: 


where  s  is  a  non-negative  scaling  coefficient;  u(t)  is  a  source  of 
amplitude  variability  of  the  mean  F;  <5(t)  is  the  random  error. 
The  function  u(t)  generates  time  variability  on  F  to  shift  the 
location  of  important  features  of  the  output  (like  the  phy¬ 
toplankton  blooms).  Representation  of  multi-variable  model 
output  (yi,y2,  ■  ■  .,ym)  can  be  found  in  Gervini  and  Gasser  (2004). 

The  problem  is  to  estimate  the  structural  mean  F.  Clearly, 
the  simple  averaging  y(t)  over  all  ensemble  realizations  under¬ 
estimates  the  amplitude  of  local  extreme  values  since  the 
peaks  vary  from  one  realization  to  another  not  only  in  inten¬ 
sity  but  also  in  timing.  To  account  for  how  it  affects  the 
structural  mean  (F),  Gervini  and  Gasser  (2004)  suggest  to 
choose  u(t)  as 

u(t)  =  u;-1(t),  (3.2) 

where  ui  is  called  the  warping  function,  which  is  represented 
by 


q 

w(t)=t+X]d^jW-  (3-3) 

1=1 

Here,  d  =  (di, . . .,  dq)  is  the  score  vector;  { ^ (t) }  are  functions 
constructed  from  a  combination  of  B-splines  and  weights  esti¬ 
mated  relative  to  a  single  model  output. 


r  =  inf(t|  z  -  z  >  e). 
t>o  1  1 


(3.5a) 


Eq.  (3.5a)  can  be  re-written  using  a  non-dimensional  toler¬ 
ance  level  e 


x  =  inf 
t>o 


(■ 


]Z-Z|  _\ 

- = -  >  E  . 

|z|  ) 


(3.5b) 


The  IPT  is  a  priori  random  value  and  depends  functionally 
on  the  reference  solution,  as  well  as  on  perturbations  in  the 
initial  conditions  and  model  parameters. 


4.  Generation  of  ensemble  perturbations 

Accuracy  in  estimating  the  non-linear  sensitivity  indices 
depends  on  ensemble  size  and  structure  of  model  paramet¬ 
rical  space.  Therefore,  to  generate  perturbations  of  model 
parameters,  it  is  better  to  get  homogeneous  coverage  in  the 
model  parameter  space.  On  the  other  hand,  for  any  model, 
the  ensemble  size  is  limited  due  to  computer  capability. 

The  Latin  hypercube  (LHC)  design  strategy  (Latin 
Hypercube,  2001)  is  used  to  generate  appropriate  pertur¬ 
bations.  An  extensive  review  of  Latin  hypercube  sampling 
technique  can  be  found  in  Helton  and  Davies  (2003);  also 
see  Rose  and  Smith  (1998)  as  an  example  of  the  ecological 


ECOLOGICAL  MODELLING  206  (2007)369-382 


373 


applications  of  this  technique.  LHC  design  increases  con¬ 
siderably  the  accuracy  in  estimating  probability  density 
functions  in  comparison  with  the  classical  Monte-Carlo 
method  (Dowling  et  al.,  1985).  For  the  same  homogeneity 
of  coverage  in  model  parametrical  space,  the  Monte-Carlo 
method  needs  an  ensemble  with  nk  components,  while  the 
LHS  strategy  only  needs  2 n(k  + 1)  components.  Here,  n  is  the 
number  of  perturbations  for  a  parameter  and  k  is  the  number 
of  model  parameters. 

Another  computational  problem  is  to  calculate  Uj  and  U_j. 
Clearly,  Eqs.  (2.5)  and  (2.10)  are  computationally  impractical.  In 
a  Monte-Carlo  frame,  it  implies  a  double  loop:  the  inner  one  is 
to  compute  E2  and  the  outer  one  is  to  compute  the  integral  over 
da,  .  Saltelli  et  al.  (1993)  suggest  two  model  parameter  matrices 
Mi  and  M2, 


(  an 

an  . 

■  alk  \ 

/  a'u 

a'12  ■ 

■  Qife\ 

Mi  = 

021 

022  ■ 

a2k 

,  m2  = 

a21 

fl22  ■ 

■  a2k 

\  0-nl 

Qn2  ■ 

■  ank  / 

\Qnl 

a'n2  ■ 

■  ank/ 

(4.1) 


Two  perturbed  parameter  matrices  are  used  for  the  specific 
parameter  oj .  The  first  one  is  constructed  from  M2  by  replacing 
the  column  [al]  to  [arj]  from  Mi,  and  another  one  from  Mi  by 
replacing  the  column  [ar)]  to  [al]  from  M2: 


(  a'll 

°12 

■  alj  ■  ■ 

■  aik\ 

a21 

a22 

■  a2j  ■  ■ 

■  a2k 

W 

<2  ■ 

anj  ■  ■ 

■  a'nkJ 

f  an 

ai2 

■  a'lj  ■ 

■  alfe  ^ 

a2i 

022 

■  «2J- 

•  a2k 

^  anl 

^n2 

.  a' .  . . 
nj 

^nk  J 

Either  Mi  or  M2  is  used  to  estimate  E(y), 

n  n 

%)  =  °r2’  ■■■■  =  a'r2’  ■■■’  ark)'  (4'3) 

r= 1  r=l 


Both  Mi  and  Nj  are  used  to  estimate  Uj, 

n 

1 


Uj  = 


■  nl>' 


>  ar(j- 1)’  arj>  ar(j+l)’ 


•  ’  ^rk) 


xf(ar  1’ - at/’  ar(j+l)’ - arfe) 


(4.4) 


Thus,  the  computational  cost  associated  with  the  full  set 
of  first-order  indices  S;  is  only  n(fe  + 1).  One  set  of  n  evaluations 
of/  is  necessary  to  compute  E,  and  k  sets  of  n  evaluations  of 
/  are  necessary  for  the  calculation  of  Uj.  Additional  k  sets  of  n 
evaluations  of  /  are  necessary  to  calculate  U_j  using  both  Mi 


and  Ny: 

1  n 

U-j  =  - - 1  y>rl,  •  ■  ■  •  ar0‘— 1) ’  arj-  ar(j+l).  ■  ■  ■  •  Qrfe) 

r—1 

x  ,  &r(j— 1),  Qfjf  ^r(j+ 1),  •  -  ■  ,  &rk)-  (4-5) 


5.  Model  of  annual  phytoplankton  cycle  for 
the  black  sea 


A  three-component  model  of  annual  phytoplankton  cycle  in 
the  Black  Sea  is  used  to  illustrate  the  non-linear  sensitivity 
analysis.  This  model  has  a  simple  configuration  and  may  not 
be  able  to  reproduce  the  phytoplankton  behavior  in  summer 
when  the  mixed  layer  depth  is  shallow.  If  the  mixed  layer 
retreats  towards  the  surface  in  the  Black  Sea,  most  of  the 
production  takes  place  below  the  mixed  layer  in  the  deep 
chlorophyll  maximum  zone. 

However,  the  model  is  able  to  reproduce  two  phyto¬ 
plankton  blooms  observed  in  reality.  The  first  one  occurs 
during  the  early  spring.  The  second  bloom  takes  place  during 
September-October.  Both  blooms  are  clearly  identified  from 
the  climatic  chlorophyll  data  (Chu  et  al.,  2005)  and  color  satel¬ 
lite  observations  (Oguz  et  al.,  2002).  This  three-component 
model  has  20  model  parameters.  The  choice  of  such  a  model  is 
not  essential  in  illustrating  the  non-linear  sensitivity  analysis. 
More  complex  biological  model  may  also  apply.  The  well- 
known  NPZ  model  originally  suggested  by  Evans  and  Parslow 
(1985)  for  the  North  Atlantic,  was  modified  and  applied  for 
reproducing  annual  phytoplankton  cycle  in  the  Black  Sea.  The 
governing  model  equations  are 


f  W .  t  (t)  =  max  [£(t) ,  0] ,  (5.1) 

dP 
dt 


a( t,  M,  P)N  c(P0  -  P)Z  m  +  f +(t) 

j  +  N  r  “  K  +  P0  -  P  M  ’ 


dz  fC(p q-p)z  m7 

dt  K  +  P0  -  P  9  M 


(5.3) 


dN 

dt 


a(t,  M,  P)N 
j  +  N 


m  +  ?+(t) ,  , 

P+ — ±A1{n0-N), 


(5.4) 


Here,  the  model  variables  are  phytoplankton  (P),  herbivore 
zooplankton  (Z),  dissolved  nutrients  (N),  all  of  them  expressed 
in  terms  of  nitrogen  specific  amount  (mmol/m3);  the  upper 
mixed  layer  depth  (M).  $(t)  is  the  time  rate  change  of  the  mixed 
layer  depth. 

The  model  parameters  are  a(t,  M,  P)  is  the  photosynthetic 
rate  of  phytoplankton,  (j ,  r)  the  half-saturation  and  mortality 
rates  for  the  phytoplankton,  (c,  K)  the  maximum  grazing  and 
half-saturation  rates  of  herbivore  zooplankton,  Po  the  phyto¬ 
plankton  threshold,  m  the  diffusivity,  /  the  grazing  efficiency, 
g  the  zooplankton  mortality  and  N0  is  the  nutrient  stock  just 
below  the  mixed  layer.  Among  these  parameters,  ^(t,  M,  P)  and 
No,  need  more  description. 

Following  Jassby  and  Platt’s  (1976)  and  Evans  and  Parslow’s 
(1985)  treatments,  the  photosynthetic  rate  of  phytoplank- 
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Table  1 

-  Model  parameters 

i 

Parameter 

Unit 

Model  Value 

1 

Nutrient  stock 

Ai 

mmol  N  m-4 

0.12 

2 

parameterization 

Ai 

mmolNm-3 

-0.75 

3 

Initial  slope  of  P/I  curve 

ao 

(Wm-2)'1  day-1 

0.025 

4 

Attenuation  coefficient  for  seawater 

kw 

m-1 

0.08 

5 

Coefficient  of  phytoplankton  self-shading 

kc 

mmolNm-3 

0.12 

6 

Half-saturation  rate  for  phytoplankton 

j 

mmolNm-3 

0.5 

7 

Mortality  rate  for  phytoplankton 

r 

day-1 

0.045 

8 

Phytoplankton  threshold 

Po 

mmolNm-3 

0.1 

9 

Grazing  efficiency  coefficient 

f 

day-1 

0.75 

10 

Zooplankton  mortality 

9 

day-1 

0.15 

11 

Phytoplankton  maximum  grow  rate 

day-1 

2 

12 

Cloudiness 

b 

— 

0.5 

Fig.  1  -  Reference  solution  and  climatic  input  for  the  Black  Sea  NPZ  model  (5.1)— (5.5):  (a)  annual  cycle  of  phytoplankton  (the 
reference  solution),  (b)  mixed  layer  depth  (given  by  Hydrometeorology  and  Hydrochemistry,  1991),  (c)  nutrient  stock  below 
the  mixed  layer  and  (d)  the  solar  irradiation  (dashed  curve)  and  photosynthetically  active  radiation  (dotted  curve). 


ton,  ce(t,  M,  P)  incorporates  several  additional  parameters, 
namely,  cloudiness  (b),  light  attenuation  by  seawater  (few)  and 
phytoplankton  self-shading  coefficient  (fec),  and  maximum 
photosynthetic  rate  (Q). 

The  nutrient  stock  (N0)  is  parameterized  by 

N0=Aiz  +  A2,  (5.5) 

which  accounts  for  vertical  variability  of  the  nutrient  distribu¬ 
tion  that  coincides  with  the  Black  Sea  climatology.  Here,  z  is 
the  vertical  coordinate. 

The  model  parameters  (Table  1)  and  the  depth  of  mixed 
layer  M  (Fig.  la)  are  taken  either  from  climatological  data 
(Hydrometeorology  and  Hydrochemistry,  1991)  or  as  tradi¬ 
tional  values  generally  accepted  in  marine  biological  modeling 
(Oguz  et  al.,  1996;  Kantha,  2004).  Fig.  lb  and  c  show  the  nutrient 
stock  below  the  mixed  layer  and  the  solar  irradiation  calcu¬ 
lated  from  the  climatological  values  of  M,  respectively. 


For  the  chosen  model  parameters  and  climatological  input 
(Fig.  la-c)  the  NPZ  model  (5.1)— (5.4)  predicts  the  existence  of 
a  bimodal  structure  with  occurrence  of  two  phytoplankton 
blooms  in  April  and  September  (Fig.  Id),  that  agrees  well  with 
the  chlorophyll-a  behavior  in  the  Black  Sea  (Chu  et  al.,  2005). 
Such  a  seasonally  varying  solution  is  treated  as  the  reference 
solution  for  the  sensitivity  studies. 


6.  Experimental  design 

During  numerical  experiments,  the  NPZ  model  (5.1)-(5.4)  is 
integrated  with  the  set  of  twelve  perturbed  parameters  (±50% 
from  the  climatic  values,  see  Table  2).  All  the  parameters  are 
uniformly  distributed  within  the  given  ranges.  Homogeneous 
coverage  in  the  model  parameter  space  for  matrices  Mi,  M2, 
Nj,  N_j  [(4.1)  and  (4.2)]  is  required  to  generate  statistically  sig¬ 
nificant  ensemble.  Let  n  (usually  nf»103)  be  the  number  of 
integrations  in  a  statistical  ensemble  for  one  parameter.  Then, 
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Table  2 

-  Ranges  of  parameter  variation 

i 

Parameter 

Range  of  Variation 

1 

Nutrient  stock 

Ai 

0.06-0.18 

2 

parameterization 

a2 

-0.375  to  -1.125 

3 

Initial  slope  of  P/I  curve 

ao 

0.01-0.04 

4 

Attenuation  coefficient  for  seawater 

few 

0.04-0.12 

5 

Coefficient  of  phytoplankton  self-shading 

fee 

0.06-0.18 

6 

Half-saturation  rate  for  phytoplankton 

j 

0.25-0.75 

7 

Mortality  rate  for  phytoplankton 

r 

0.0225-0.0675 

8 

Phytoplankton  threshold 

Po 

0.05-0.15 

9 

Grazing  efficiency  coefficient 

/ 

0. 5-1.0 

10 

Zooplankton  mortality 

g 

0.075-0.225 

11 

Phytoplankton  maximum  grow  rate 

1-3 

12 

Cloudiness 

b 

0.25-0.75 

for  the  NPZ  model  with  k  =  12,  only  26,000  perturbations  are 
generated.  This  is  computationally  feasible  even  for  a  personal 
computer. 

There  are  about  20  parameters  in  original  Evans  and 
Parslow  (1985)  biological  model,  but  only  12  parameters  here 
for  the  sensitivity  study.  From  the  computational  point  of  view, 
there  is  no  restriction  on  the  number  of  perturbed  parame¬ 
ters.  However,  part  of  model  parameters  are  fixed  due  to:  (a) 
little  influence  on  the  model  output  (non-control  parameters) 
and  (b)  little  information  on  the  range  of  variation  for  that 
parameter. 


7.  Construction  of  warping  functions 

Since  the  reference  solution  is  quasi-periodic,  one  year  is 
taken  for  constructing  the  warping  functions  for  annual 
phytoplankton  cycle.  Our  computations  show  that  a  two- 
component  decomposition  [i.e.,  q  =  2  in  (3.3)]  is  sufficient  for 
the  reconstruction  of  the  mean  phytoplankton  annual  cycle 


shown  in  Fig.  Id  from  the  ensemble  of  perturbed  reference 
solutions  (Fig.  2a).  This  ensemble  of  stochastic  realizations 
is  calculated  by  adding  white  noise  to  the  reference  solution 
(Fig.  Id).  Here,  ratio  of  noise-to-signal  (the  reference  solu¬ 
tion)  is  about  1.  Two  bell-shaped  warping  functions  i^i  and 
(Fig.  2b)  are  computed  on  the  base  of  six  B-splines  with  equally 
spaced  knots. 

Reconstruction,  using  fi  and  f2,  leads  to  the  annual 
phytoplankton  cycle  with  the  maxima  corresponding  approx¬ 
imately  to  the  time  of  the  two  blooms  (Fig.  2c).  Some 
distortions  of  the  reconstructed  reference  solution  (1-month 
shift  of  the  fall  bloom  and  reduction  of  bloom  peaks) 
observed  in  Fig.  2c  are  due  to  high  noise  added  to  the 
reference  solution.  However,  the  accuracy  of  the  reconstruc¬ 
tion  increases  with  reducing  the  noise-to-signal  ratio.  The 
root  mean  square  error  indicates  that  this  technique  can 
reconstruct  the  two  modal  structures  in  the  annual  phy¬ 
toplankton  cycle  with  reasonable  accuracy  using  only  two 
warping  functions  even  with  large  noise-to-signal  ratios 
(Fig.  2d). 


Fig.  2  -  Construction  and  evaluation  of  warping  functions:  (a)  an  ensemble  of  perturbed  reference  solutions,  (b)  warping 
functions  and  1//2,  (c)  reconstructed  mean  and  (d)  the  root  mean  square  error  between  the  reference  solution  and  the 
reconstructed  mean  solution  f  2  estimated  by  the  warping  function  with  q  =  2. 


376 


ECOLOGICAL  MODELLING  206  (2OO7)  369-382 


Fig.  3  -  The  first-order  sensitivity  of  phytoplankton  amount 
in  upper  mixed  layer  to  variations  of  twelve  model 
parameters.  Labels  1, . . .,  12  correspond  to  parameters 
specified  as  in  Tables  1  and  2. 


Fig.  4  -  The  high  order  sensitivity  of  phytoplankton  amount 
in  upper  mixed  layer  to  variations  of  twelve  model 
parameters.  Labels  1,  . . .,  12  correspond  to  parameters 
specified  as  in  Tables  1  and  2. 


8.  Sensitivity  indices  for  phytoplankton 
component 

The  above  non-linear  sensitivity  analysis  is  directly  applied 
to  phytoplankton  concentration  computed  by  the  NPZ  model. 
The  first-order  sensitivity  indices  (Sj)  for  all  12  parameters 
(Table  2)  are  calculated  numerically  for  each  month  within 
a  year  cycle.  Fig.  3  shows  the  temporal  variations  of  Sj  = 
Y'_, S,-(i  =  1 _ 12)  (indicated  by  bold  curves  in  Fig.  3). 

Two  interesting  features  are  important  for  further  analy¬ 
sis.  First,  the  set  of  influential  parameters  considerably  varies 
with  time.  The  annual  phytoplankton  cycle  is  divided  into 
three  periods:  (1)  prior  to  spring  bloom,  (2)  between  the  spring 
and  fall  blooms  and  (3)  after  the  fall  bloom.  It  is  expected  that 
the  phytoplankton  during  the  first  and  third  periods  is  mostly 
influenced  by  initial  slope  of  the  P/I  curve  a0  (i  =  3)  and  the 
light  conditions  (i  =  4)  while  the  phytoplankton  during  the  sec¬ 
ond  period  is  driven  mainly  by  the  nutrients  stock  under  the 
mixed  layer  (i  =  1,  2). 

Second,  the  summation  of  the  first-order  indices  never 
reaches  1.  This  means  that  the  model  output  is  not  simple 
additive  relative  to  parameters.  Even  if  a  possible  uncertainty 
up  to  50%  in  the  model  output  is  taken  into  account,  the  sum 
of  all  first  indices  drops  evidently  during  the  blooms.  It  means 
that  the  contribution  of  interactions  among  model-related 
errors  generated  by  perturbations  of  different  parameters,  to 
model  output  is  large.  The  first-order  indices  alone  cannot 
fully  explain  the  model  output  variance  as  a  whole. 

Contribution  of  interactions  between  model-related  errors 
generated  by  perturbations  of  different  parameters  is  esti¬ 
mated  by  the  difference, 


Fig.  4  shows  the  temporal  variation  of  Sjnt  = 

Although,  we  can  still  conclude  higher  significance  of  some 
factors,  such  as  the  initial  slope  a  of  the  I/P  curve  and  the 
parameters  defining  the  nutrient  stock  (Ai,  A2)  below  the 
mixed  layer,  the  interpretation  of  this  plot  is  questionable. 
Therefore,  instead  of  analyzing  the  temporally  varying  model 


output,  we  use  the  scale  parameters  and  scores  of  the  corre¬ 
sponding  functions  as  described  above. 


9.  Sensitivity  indices  for  scale  (s)  and  score 
vector  (di,  d2) 

Section  7  shows  that  the  model  output  (phytoplankton  con¬ 
centration)  can  be  represented  only  by  two  warping  functions 
(fi,  i'l)-  Therefore,  the  model  output  sensitivity  is  reduced  to 
the  sensitivity  of  a  scale  s  and  a  score  vector  (di,  d2)  relative  to 
variations  of  model  parameters. 

9.1.  Sensitivity  0/ di 

The  first  component  of  score  vector  (di)  in  general,  determines 
the  spring-bloom  landmark  and  fixes  the  position  of  the  spring 
bloom.  Therefore,  joint  analysis  of  the  single  (Fig.  5a)  and 
total  (Fig.  5b)  sensitivity  indices  for  di  identifies  the  main  fac¬ 
tors  responsible  for  the  generation  of  this  bloom.  The  indices 
(Fig.  5a  and  b)  are  computed  with  95%  confidence  interval 
using  the  bootstrap  technique  (Efron  and  Tibshirani,  1993) 
with  a  bootstrap  sample  dimension  of  103. 

It  is  obvious  from  Fig.  5a  and  b  that  the  single  indices 
are  very  small  as  comparing  to  the  total  indices.  This 
implies  high  importance  of  parameter  interactions.  The  sin¬ 
gle  indices  for  di  show  that  the  initial  slope  a0  of  the  P/I 
curve  (i  =  3)  has  a  significant  effect  (about  0.12)  while  all  the 
other  parameters  are  non-significant  (Fig.  5a).  This  result 
seems  to  be  correct  only  for  small  perturbations  of  model 
parameters. 

In  contrast  to  this,  the  total  indices  for  di  give  more  infor¬ 
mation  about  model  sensitivity  and  factors  determining  it. 
Fig.  5b  indicates  that  the  parameter  a 0  (i  =  3)  still  retains  the 
leading  role  with  the  total  index  of  about  1,  however,  other 
parameters,  such  as  the  light  attenuation  coefficient  of  the 
seawater  (i  =  4)  and  cloudiness  (i  =  12)  become  important  too. 
The  total  indices  for  these  parameters  are  up  to  0.85-0.87.  The 
light  attenuation  coefficient  and  cloudiness  are  related  mostly 
to  the  light  limitation  of  the  phytoplankton  growth. 
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Fig.  5  -  (a)  Single  and  (b)  total  sensitivity  indices  for  the  first  score  vector  component  (dj).  Bars  show  95%  confidence  interval 
computed  by  the  bootstrap  technique. 


Less  important  but  still  influential  is  the  mortality  rate 
of  the  phytoplankton  (i  =  7)  and  the  parameter  Ai  determin¬ 
ing  the  nutrient  stock  below  the  mixed  layer  (i  =  1).  The  total 
indices  for  these  parameters  reach  0.61-0.65.  The  model  solu¬ 
tion  is  not  sensitive  to  the  choice  of  second  parameter  A2  in 
Eq.  (5.4)  (i  =  2),  half-saturation  rate  for  phytoplankton  (i  =  6) , 
crazing  efficiency  coefficient  (i  =  9)  and  zooplankton  mortal¬ 
ity  (i  =  10).  The  weakest  sensitivity  is  also  for  the  coefficient 
of  plankton  self-shading  (i  =  5),  phytoplankton  threshold  (i  =  8) 
and  phytoplankton  maximum  growth  rate  (i  =  11). 

Our  method  estimates  relative  importance  of  each  model 
parameters  on  the  spring  bloom  generation  and  identi¬ 
fies  their  contributions  using  non-dimensional  numbers  j/j  = 
ST /maxfS^)  as  follows 
i 

0.85  <  (i)3,  i)4,  i)12)  <  1.00  for  highly  significant, 

0.60  <  (i/i,  1/7)  <  0.65  for  significant, 

0.38  <  (i)2,  i)6,  i/9,  r)  10)  <  0.45  for  low  significant, 

0.15  <  (1)5,  i)g,  »)n)  <  0.25  for  not  significant. 

Sensitivity  enhancement  of  the  spring  bloom  (phy¬ 
toplankton  growth)  to  the  light  limitation  is  intuitively 
understandable.  However,  ranging  model  parameters  along 
their  contributions  to  model  output  variability  does  not  follow 
exactly  the  physical  intuition. 

More  biological  insight  can  be  obtained  using  the  second- 
order  indices  Sjj  (Fig.  6).  These  indices  show  that  the  sensitivity 
of  model  output  depends  strongly  on  interactions  among 
model-related  errors  caused  by  different  model  parameters. 
For  example,  the  model-related  error  due  to  uncertainty  of  ini¬ 
tial  slope  of  the  P/I  curve  (i  =  3)  is  amplified  through  interaction 
with  error  caused  by  uncertainty  of  phytoplankton  threshold 
(i  =  8)  (large  value  of  S38  in  Fig.  6),  and  is  diminished  due  to 
uncertainty  of  grazing  efficiency  coefficient  (i  =  9)  (small  value 
of  S39  in  Fig.  6).  However,  the  model  output  is  very  sensitive 
to  uncertainty  in  determination  of  parameter  pairs  (i,j)  =  (2,3), 
(2,10),  (4,6),  (5,10)  and  (11,12).  In  general,  we  may  conclude 
that  the  generation  of  the  spring  bloom  depends  strongly  on 
light-related  model  parameters  and  weakly  on  different  phy¬ 
toplankton  ratios. 


9.2.  Sensitivity  0/ d2 

The  second  component  of  score  vector,  d2,  determines  the 
fall-bloom  landmark.  The  single  sensitivity  indices  (Fig.  7a) 
of  d2  show  that  the  nutrient  stock  parameterization  (i  =  1,  2) 
is  the  most  important  for  reproducing  the  fall  phytoplank¬ 
ton  bloom  (Si  =  0.12  and  S2  =  0.10).  On  the  contrary,  the  total 
sensitivity  indices  identify  six  most  important  model  param¬ 
eters  for  reproducing  the  fall  phytoplankton  bloom  (Fig.  7b). 
They  are  the  nutrient  stock  parameterization  (i  =  1),  the  ini¬ 
tial  slope  of  the  P/I  curve  (i  =  3) ,  the  attenuation  coefficient  for 
seawater  (i  =  4),  the  mortality  rate  for  phytoplankton  (i  =  7),  the 
zooplankton  mortality  (i  =  10)  and  the  cloudiness  (i  =  12).  The 
reference  solution,  in  general,  is  insensitive  to  perturbations 
of  parameters  with  i  =  2,  5,  6,  8,  9  and  11.  Significant  total  effect 
for  the  fall-bloom  landmark  also  depends  on  the  parameters 
for  the  nutrient  distribution  below  the  mixed  layer  as  well  as 
photosynthesis. 

The  sensitivity  indices  Sjj  (Fig.  8)  show  the  important 
effect  of  interactions  among  model-related  errors  generated 
by  uncertainty  inserted  into  the  following  pairs  of  model 
parameters  (ij)  =  (l,3),  (1,4),  (1,12),  (1,  6-11),  (3,4),  (3,11),  (6,9), 
(7,10),  (8,9)  and  (9,11).  Again,  note  that  interactions  among 
model  related  errors  can  lead  to  increasing  or  decreasing 


Fig.  6  -  Second-order  indices  Sy  for  the  first  component  of 
the  score  vector  di. 
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Fig.  7  -  (a)  Single  and  (b)  total  sensitivity  indices  for  the 
second  score  vector  component  (d2).  Bars  show  95% 
confidence  interval  computed  by  the  bootstrap  technique. 


model  output  sensitivity  to  uncertain  model  parameters. 
For  example,  the  reference  solution  is  insensitive  to  per¬ 
turbations  inserted  into  model  parameters  pairs  (i,j)  =  (l,7), 
(1,8),  (1,9),  (1,10),  (1,11).  The  large  value  of  single  sensitivity 
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Fig.  8  -  Second-order  indices  Si;-  for  the  second  component 
of  the  score  vector  d2. 
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index  for  Aj  indicates  high  influence  of  Ai  on  the  reference 
solution. 

9.3.  Sensitivity  of  the  scale  coefficient  s 

The  single  sensitivity  indices  of  the  scale  coefficient  (s)  demon¬ 
strate  no  one  parameter  (or  several  parameters)  dominating 
the  model  output  (Fig.  9a).  Although  Si  and  S2  are  larger 
than  other  indices,  their  values  are  quite  small  (<0.11).  How¬ 
ever,  the  total  sensitivity  indices  (Fig.  9b)  and  second-order 
indices  Si;  (Fig.  10)  indicate  high  sensitivity  of  the  scale  coeffi¬ 
cient  (s)  to  interactions  among  model-related  errors  caused  by 
uncertainty  inserted  in  different  pairs  of  model  parameters. 
The  model  parameters  are  decomposed  into  two  groups  with 
respect  to  model  output  sensitivity.  The  first  group  contains 
model  parameters  relative  to  which  a  has  maximum  sensi¬ 
tivity:  nutrient  stock  parameterization  (i  =  1,  2),  initial  slope  of 
P/I  curve  (i  =  3)  and  attenuation  coefficient  for  seawater  (i  =  4). 
The  second  group  includes  model  parameters  causing  mini¬ 
mum  sensitivity  of  model  output:  from  plankton  self-shading 
(i  =  5)  to  cloudiness  (i  =  12).  The  second-order  indices  S;j  (Fig.  10) 
show  that  the  model  output  is  strongly  sensitive  to  interac¬ 
tions  among  model-related  errors  induced  by  pairs  of  model 
parameters:  (1,2),  (1,3),  (1,4),  (1,10)  and  (2,3). 


Model  Parameters 


Fig.  9  -  (a)  Single  and  (b)  total  sensitivity  indices  for  the  scaling  coefficient  s.  Bars  show  95%  confidence  interval  computed 
by  the  bootstrap  technique. 


ECOLOGICAL  MODELLING  206  (2007)369-382 


379 


-  8 


10.07 


0.06 


0.04 


0.03 


0  01 


1  23  456  7  89  10  11 


Fig.  10  -  Second-order  indices  Sy  for  the  scaling  coefficient  s. 


10.  Sensitivity  Indices  of  IPT 

The  warping  functions  are  good  metrics  for  the  non-linear 
sensitivity  analysis  if  the  model  output  can  be  represented  by 
few  warping  functions.  This  only  happens  when  the  reference 
solution  presents  the  periodical  or  quasi-periodical  attractor. 
For  more  complex  dynamical  regimes,  such  as  chaotic  attrac¬ 
tor,  large  number  of  the  warping  functions  are  needed  for 
decomposition  of  model  output.  To  overcome  this  difficulty, 
IPT  (represented  by  r)  is  used  to  quantify  model  non-linear 
sensitivity.  Clearly,  smaller  values  of  IPT  correspond  to  higher 
sensitivity  and  vice  versa. 

Analyzing  the  functional  dependence,  r  =  z(s),  we  can 
understand  the  sensitivity  of  the  reference  solution  to  pertur¬ 
bations  of  various  intensities.  The  tolerance  level  s  controls 
intensities  of  model-related  errors  caused  by  perturbations  in 
model  parameters. 

The  IPT,  in  general,  is  a  non-smooth  function  in  12- 
dimensional  parametrical  space  (ai, . . .,  012).  Fig.  11  shows  the 
IPT  computed  in  the  model  parameter  sub-space  [a\,  £12)  for 
two  tolerance  levels  with  high  “ridge”  and  deep  “valleys”.  IPT 
depends  strongly  not  only  on  magnitude  but  also  on  direction 
of  the  perturbation  vector. 

The  analysis  of  the  IPT  landscapes  between  £  =  0.07 
(Fig.  11a)  and  e  =  0.2  (Fig.  lib)  demonstrates  that  the  model- 
related  errors  do  not  grow  in  dependently  for  tolerance  levels 


t (days) 


Fig.  11  -  IPT  for  the  NPZ  model  (5.1)-(5.4)  in  the  model 
parameter  phase  sub-space  (Ai,  A2)  with  the  tolerance 
level:  (a)  e  =  0.07  and  (b)  e  =  0.2. 


larger  than  0.15-0.20.  It  indicates  that  the  errors  interact 
among  themselves,  and  the  total  effect  from  the  parameter 
perturbations  on  model  output  is  non-linear.  For  small  tol¬ 
erance  levels  (as  an  example,  £  =  0.07),  the  single  indices  (Sj) 
are  small  (Fig.  12a).  Maximum  value  for  the  indices  is  0.15.  It 
indicates  weak  sensitivity  of  the  model  output  to  the  choice 
of  model  parameters.  However,  the  total  indices  (Sj)  show 


Fig.  12  -  (a)  Single  and  (b)  total  sensitivity  indices  for  the  IPT  with  the  tolerance  level  of  e  =  0.07.  Bars  show  95%  confidence 
intervals  computed  by  the  bootstrap  technique. 
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Fig.  13  -  Second-order  indices  Sjj  for  the  IPT  with  the 
tolerance  level  of  e  =  0.07. 


at  least  seven  influential  parameters  (Fig.  12b):  the  nutrient 
stock  parameterization  (i  =  1),  the  initial  slope  of  P/I  curve  (i  =  3), 
the  attenuation  coefficient  for  seawater  (i  =  4),  the  coefficient 
of  phytoplankton  self-shading  (i  =  5),  the  half-saturation  rate 
for  phytoplankton  (i  =  6),  the  mortality  rate  for  phytoplankton 
(i  =  7)  and  the  cloudiness  (i  =  12).  Obviously,  the  biological  NPZ 
model  (5.1)-(5.4)  is  not  good  enough  for  long-term  prediction 
with  high  accuracy  because  it  is  difficult  to  identify  the  seven 
parameters  from  data  with  appropriate  accuracy.  The  model 
produces  the  right  type  of  attractor  but  cannot  reproduce  such 
fine  details  as  the  phase  of  oscillations.  Note  that  interactions 
among  the  model-related  errors  caused  by  different  param¬ 
eter  perturbations  do  not  add  to  the  general  prediction  error 
(Fig.  13).  In  this  regime,  the  model  sensitivity  can  be  analyzed 
by  a  tangent  model,  i.e.,  through  the  first-order  sensitivity 
functions. 

For  intermediate  values  of  i  ranging  from  0.1  to  0.5,  the 
total  sensitivity  indices  identify  the  model  parameters  with 
maximum  contribution  to  model  output  sensitivity  (Fig.  14a 
and  b).  The  reference  solutions  are  most  sensitive  to  the  choice 
of  the  initial  slope  of  P/I  curve  (i  =  3),  the  attenuation  coefficient 
for  seawater  (i  =  4)  and  the  cloudiness  (i  =  12) . 

The  sensitivity  of  model  output  to  variations  of  the  nutrient 
stock  parameterization  (i  =  1,  2)  and  different  phytoplank¬ 
ton  ratios  is  weak.  However,  the  second-order  indices  Sj j 
(Fig.  15)  show  that  the  contribution  of  the  nutrient  stock 
parameterization  to  the  model  output  sensitivity  can  be 
amplified  considerably  through  interactions  among  appro¬ 
priate  model-related  errors.  For  example,  the  model-related 
errors  corresponding  to  the  perturbations  of  Ai  and  A2 
strongly  interact  with  the  errors  caused  by  perturbations  in 
ao  and  kw.  These  interactions  are  non-linear  and  do  not  exist 
for  small  tolerance  levels  (see  Figs.  12b  and  13). 

Thus,  a  model  of  the  Black  Sea  phytoplankton  cycle  with 
higher  accuracy  than  the  NPZ  model,  requires  high  accurate 
model  parameters.  It  is  hard  to  estimate  these  parameters 
with  necessary  accuracy  from  existing  biological  observations. 
A  less  accurate  model,  such  as  the  NPZ  model,  which  is  able 
to  reproduce  the  main  events  (spring  and  fall  blooms),  may  be 
successfully  constructed  because  only  few  parameters  need 
to  be  identified. 


Fig.  14  -  The  total  sensitivity  indices  for  the  IPT  on  the 
tolerance  level:  (a)  i  =  0.2  and  (b)  i  =  0.5.  Bars  show  95% 
confidence  intervals  computed  by  the  bootstrap  technique. 


i 

Fig.  15  -  Second-order  indices  S y  for  the  IPT  with  the 
tolerance  level  of  e  =  0.2. 


11.  Conclusions 

The  present  study  has  developed  a  variance-based  method 
for  the  analysis  of  non-linear  sensitivity  of  marine  biological 
models  to  large  variations  of  model  parameters.  This  approach 
utilizes  the  Sobol’-Saltelli’s  sensitivity  indices  as  a  measure 
of  model  sensitivity  and  special  scalar  model  output  decom- 
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posed  through  the  warping  functions  or  represented  by  the 
irreversible  predictability  time. 

The  method  does  not  identify  directly  the  optimal  param¬ 
eters  from  model-data  fitting  but  ranges  the  parameters 
relative  to  their  influence  to  model  output.  It  determines  the 
control  parameters,  which  contribute  most  to  output  variabil¬ 
ity.  Additional  research  is  required  to  increase  knowledge  of 
the  control  parameter  in  order  to  reduce  output  uncertainty. 
The  insignificant  parameters  are  held  constant  or  even  elimi¬ 
nated  from  the  model. 

Many  methods  currently  in  use  have  some  sort  of  sensi¬ 
tivity  aspect  for  linear  and  non-linear  biological  models.  The 
variance-based  method  developed  here,  can  be  used  to  investi¬ 
gate  non-linear  sensitivity  of  model  output  to  large  variations 
of  model  parameters. 

Large  values  of  the  sensitivity  indices  S;,-  show  the  amplifi¬ 
cation  of  the  model  output  variance  due  to  perturbations  of 
both  ith  and  jth  model  parameters.  Here,  output  variances 
corresponding  to  perturbations  of  appropriate  parameters  are 
not  additive.  Contrarily,  small  values  of  the  indices  S;;  show 
no  such  amplification  of  the  model  output  variance  when  two 
parameters  are  perturbed  simultaneously. 

The  simplified  three-component  model  of  annual  phyto¬ 
plankton  cycle  in  the  Black  Sea  is  chosen  to  illustrate  the 
technique  and  to  understand  some  generic  features  of  model 
sensitivity  to  large  perturbations  of  model  parameters.  This 
model  has  very  simple  configuration  and  may  fail  to  repro¬ 
duce  the  phytoplankton  behavior  in  summer  when  the  mixed 
layer  depth  is  shallow.  However,  the  model  reproduces  the  bi- 
modal  behavior  of  phytoplankton  observed  in  the  Black  Sea 
and  can  be  a  useful  tool  for  the  demonstration  of  capability  of 
the  non-linear  sensitivity  analysis. 

The  single  and  total  sensitivity  indices  demonstrate  that 
the  model  predicted  spring  bloom  (in  April-May)  is  most 
sensitive  to  the  choice  of  the  initial  slope  of  the  photosyn¬ 
thesis/irradiation  curve.  As  far  as  the  parameter  interactions 
are  concerned,  the  light  attenuation  of  the  seawater,  cloudi¬ 
ness,  the  mortality  rate  of  the  phytoplankton  and  the  nutrient 
stock  parameter  become  influential  as  compared  to  the  initial 
slope  of  the  P/I  curve,  the  light  attenuation  coefficient  of  the 
seawater  and  cloudiness.  These  factors  are  mostly  related  to 
the  light  limitation  of  the  phytoplankton  growth.  Of  course, 
this  result  is  intuitive  from  the  physical  point  of  view. 

However,  ranging  model  parameters  along  their  con¬ 
tributions  to  model  output  variability  does  not  show 
straightforward  result.  The  intuition  is  also  useless  when  we 
try  to  understand  and  to  estimate  contributions  of  interactions 
among  model-related  errors  caused  by  perturbations  of  differ¬ 
ent  parameters  to  model  output  variance.  For  example,  this  is 
not  a  priori  clear  when  model-related  error  caused  by  uncer¬ 
tainty  in  the  nutrient  stock  parameterization  is  amplified 
considerably  by  interactions  with  uncertainty  in  the  light- 
related  factors. 
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